## Prioritization preferences for COVID-19  -------------------------
## vaccination are consistent across five countries -----------------
## Reproduction materials (Main text) -------------------------------

set.seed(1234)

# load packages
source("code/helpers.R")

# load data
tab_2_df <- readRDS("data/tab_2_df.rds")
ranking_df <- readRDS("data/ranking_df.rds")
conjoint_df <- readRDS("data/conjoint_df.rds")
plot_helper_df <- readRDS("data/plot_helper_df.rds")

## -----------------------------------------------------------------
## Table 2. Overview of studies ------------------------------------
## -----------------------------------------------------------------

tab_2_df %>%
  gt::gt() %>%
  # Add flag images
  gtExtras::gt_img_multi_rows(columns = Country, height = 19) %>%
  gtExtras::gt_theme_nytimes() %>%
  gt::tab_row_group(rows = c(Study == "Our study"),
                    label = " ") %>%
  gt::tab_row_group(rows = c(Study != "Our study"),
                    label = "") %>%
  gt::fmt_markdown(columns = "Preference elicitation device") %>%
  gt::cols_width(Study ~ px(100),
                 Country ~ px(100),
                 `Prioritization groups/characteristics` ~ px(400),
                 `Highest priorities` ~ px(200),
                 `Preference elicitation device` ~ px(80),
                 everything() ~ px(80)) %>%
  # General table options
  tab_options(
    data_row.padding = px(2),
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.style = "hidden",
    table_body.border.top.style = "solid",
    column_labels.border.bottom.style = "solid"
  ) %>%
  tab_style(
    style = list(cell_text(color = "#000", align = "center",v_align = "middle", weight = "bold")),
    locations = cells_column_labels(columns = c(everything())
    )) %>%
  cols_align(
    align = "center",
    columns = c(`Preference elicitation device`,`Field time`)
  ) %>%
  gt::gtsave("figures/table-2-evidence-classification.png")


## -----------------------------------------------------------------
## Figure 2. Average ranks of groups -------------------------------
## -----------------------------------------------------------------

mean_rank_df <- ranking_df %>%
  group_by(country, type) %>%
  summarise(mean_rank = mean(val, na.rm=T))

order_rank <- rev(mean_rank_df %>%
                    dplyr::arrange(country, mean_rank) %>%
                    dplyr::group_by(country) %>%
                    dplyr::mutate(rank = dplyr::row_number()) %>%
                    dplyr::pull(type) %>% .[1:8])

fig_2 <- mean_rank_df %>%
  dplyr::arrange(country, mean_rank) %>%
  dplyr::group_by(country) %>%
  dplyr::mutate(rank = dplyr::row_number()) %>%
  ggplot(aes(y = mean_rank - 1, x = factor(type, levels = order_rank))) +
  geom_bar(stat = "identity") +
  geom_label(aes(y = mean_rank - 0.5, x = type, label = rank), color = "darkblue") +
  geom_text(aes(y = mean_rank - 1.75, x = type, label = paste0("(",sprintf('%.1f', mean_rank),")")), color = "#ffffff", hjust = 0, size = 2.7, fontface = "bold") +
  coord_flip() +
  facet_wrap(~country) +
  theme_bw() +
  scale_x_discrete(labels = order_rank) +
  scale_y_continuous(limits = c(-0.25,8.5),
                     labels = function(y) y + 1,
                     breaks = c(1,3,5,7)) +
  theme(legend.position = "none", text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        )) +
  labs(x = "", y = "Mean Rank") 

fig_2

ggsave("figures/fig-2-ranking-bar.png", plot=fig_2, height=9, width=12)

## -----------------------------------------------------------------
## Figure 3. Conjoint marginal means -------------------------------
## -----------------------------------------------------------------

# estimation function
est_mm <- function(df){
  # linear model
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition + Citizenship, data = df, weights = weight, clusters = personid)
  
  # gather marginal means
  vars <- c("Age", "EarlyRegistration", "Gender" ,"Children" , "Job" , "Precondition", "Citizenship")
  df_res <- do.call( rbind, lapply(vars, function(v) {  
    df_mm <-  data.frame(emmeans(m,v))
    df_mm$names <- paste0(names(df_mm)[1],df_mm[[1]])
    names(df_mm)[1] <- "var"  
    return(df_mm)
  }))
} 

# country estimates
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_mm(.))

# pooled estimate
df_res_pooled <- est_mm(conjoint_df)
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- dplyr::bind_rows(df_res_pooled, df_res)
fig_3_df <- dplyr::left_join(plot_helper_df, df_res_all) %>%
  dplyr::mutate(country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))

# ticks and labs
xaxis_ticks <- unique(fig_3_df$position)
xaxis_lab <- unique(paste0(fig_3_df$attr, " | ", fig_3_df$levl)) %>% stringr::str_replace(., "Children", "Has children")
vertical_lines <- c(5,8,11,14,17,24)

# create figure
fig_3 <- ggplot(fig_3_df, aes(x=position,y=emmean,ymin=lower.CL,ymax=upper.CL)) + 
  geom_point(size=1.3) + 
  geom_linerange(size=0.7) + 
  geom_hline(yintercept = 0.5,alpha=0.3,col="red") +
  facet_wrap(~ country)  + 
  theme_bw() +
  ylim(0.28,0.72) +
  scale_x_continuous(breaks=xaxis_ticks,
                     labels = xaxis_lab) +
  geom_vline(xintercept=vertical_lines, col="grey", linetype="dashed", size=1) +
  coord_flip() +  xlab("Vaccination Candidates' Attributes") + ylab("Marginal Mean of Prioritization for Vaccination") + 
  scale_color_grey() + 
  theme_bw() +
  theme(legend.position = "none", text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

fig_3

ggsave("figures/fig-3-margmeans.png", plot=fig_3, height=8, width=10)
